// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtSubdivideQuadsFilter.C
// ************************************************************************* //

#include <avtSubdivideQuadsFilter.h>

#include <vtkVisItUtility.h>
#include <vtkTriangleFilter.h>
#include <vtkCell.h>
#include <vtkCellData.h>
#include <vtkDataSet.h>
#include <vtkIdList.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkTemplateAliasMacro.h>

#include <avtCallback.h>

#include <string>

// ****************************************************************************
//  Method: avtSubdivideQuadsFilter constructor
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Tue Nov 2 06:04:08 PDT 2004
//
// ****************************************************************************

avtSubdivideQuadsFilter::avtSubdivideQuadsFilter()
{
    haveIssuedWarning = false;
}


// ****************************************************************************
//  Method: avtSubdivideQuadsFilter destructor
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Tue Nov 2 06:04:08 PDT 2004
//
//  Modifications:
//
// ****************************************************************************

avtSubdivideQuadsFilter::~avtSubdivideQuadsFilter()
{
}


// ****************************************************************************
//  Method:  avtSubdivideQuadsFilter::Create
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Tue Nov 2 06:04:08 PDT 2004
//
// ****************************************************************************

avtFilter *
avtSubdivideQuadsFilter::Create()
{
    return new avtSubdivideQuadsFilter();
}


// ****************************************************************************
//  Method:      avtSubdivideQuadsFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Tue Nov 2 06:04:08 PDT 2004
//
// ****************************************************************************

void
avtSubdivideQuadsFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const SubdivideQuadsAttributes*)a;

    // We need to specify that we want a secondary variable as soon as
    // possible.
    if (strcmp(atts.GetVariable().c_str(), "default") != 0)
    {
        SetActiveVariable(atts.GetVariable().c_str());
    }
}


// ****************************************************************************
//  Method: avtSubdivideQuadsFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtSubdivideQuadsFilter with the given
//      parameters would result in an equivalent avtSubdivideQuadsFilter.
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Tue Nov 2 06:04:08 PDT 2004
//
// ****************************************************************************

bool
avtSubdivideQuadsFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(SubdivideQuadsAttributes*)a);
}


// ****************************************************************************
//   Templatized helper methods.
//
// ****************************************************************************

template <class T> int
DetermineSubdivisions(vtkIdType id[4], int maxSubdivs, double threshold, T *p)
{
    int numSubdivs = 0;
    T v[4];
    v[0] = p[id[0]];
    v[1] = p[id[1]];
    v[2] = p[id[2]];
    v[3] = p[id[3]];
    T min = v[0];
    T max = v[0];
    min = (v[1] < min ? v[1] : min);
    max = (v[1] > max ? v[1] : max);
    min = (v[2] < min ? v[2] : min);
    max = (v[2] > max ? v[2] : max);
    min = (v[3] < min ? v[3] : min);
    max = (v[3] > max ? v[3] : max);
    T diff = max-min;
    while ((diff > threshold) && (numSubdivs < (maxSubdivs-1)))
    {
        diff /= 2.;
        numSubdivs++;
    }
    return numSubdivs;
}

int
DetermineSubdivisions(vtkDataArray *arr, vtkIdType id[4], int max, 
                      double threshold)
{
    int numSubdivs = 0;
    switch(arr->GetDataType())
    {
        vtkTemplateAliasMacro(numSubdivs = DetermineSubdivisions(id, max, threshold,
            static_cast<VTK_TT *>(arr->GetVoidPointer(0))));
    }
    return numSubdivs;
}

template <class T> inline void
CreateEdgePoint(vtkIdType id[4], T pt[3], T *pts_ptr)
{
    pt[0]  = 0.5*pts_ptr[3*id[0]+0];
    pt[0] += 0.5*pts_ptr[3*id[1]+0];
    pt[1]  = 0.5*pts_ptr[3*id[0]+1];
    pt[1] += 0.5*pts_ptr[3*id[1]+1];
    pt[2]  = 0.5*pts_ptr[3*id[0]+2];
    pt[2] += 0.5*pts_ptr[3*id[1]+2];
}

template <class T> inline void
CreateMidPoint(bool fanOut, bool iSpecial, bool jSpecial, 
    vtkIdType id[4], double w[4], double middle[3],
    T pt[3], T *pts_ptr)
{
    pt[0]  = w[0]*pts_ptr[3*id[0]+0];
    pt[0] += w[1]*pts_ptr[3*id[1]+0];
    pt[0] += w[2]*pts_ptr[3*id[2]+0];
    pt[0] += w[3]*pts_ptr[3*id[3]+0];
    pt[1]  = w[0]*pts_ptr[3*id[0]+1];
    pt[1] += w[1]*pts_ptr[3*id[1]+1];
    pt[1] += w[2]*pts_ptr[3*id[2]+1];
    pt[1] += w[3]*pts_ptr[3*id[3]+1];
    pt[2]  = w[0]*pts_ptr[3*id[0]+2];
    pt[2] += w[1]*pts_ptr[3*id[1]+2];
    pt[2] += w[2]*pts_ptr[3*id[2]+2];
    pt[2] += w[3]*pts_ptr[3*id[3]+2];
    if (fanOut && (iSpecial || jSpecial))
    {
        T diff[3];
        diff[0] = pt[0] - middle[0];
        diff[1] = pt[1] - middle[1];
        diff[2] = pt[2] - middle[2];
        pt[0] += 0.01*diff[0];
        pt[1] += 0.01*diff[1];
        pt[2] += 0.01*diff[2];
    }
} 


// ****************************************************************************
//  Method: avtSubdivideQuadsFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the SubdivideQuads filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Tue Nov 2 06:04:08 PDT 2004
//
//  Modifications:
//    Kathleen Biagas, Thu Aug  9 10:29:58 PDT 2012
//    Use new templatized methods for a few things, to handle multiple data 
//    types for vars and coords.
//
//    Eric Brugger, Tue Aug 19 09:22:22 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtSubdivideQuadsFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    //
    // Get the information about how this filter should execute.
    //
    double threshold  = atts.GetThreshold();
    int    maxSubdivs = atts.GetMaxSubdivs()+1; // +1 makes indexing easier
                                                // and allows user to have
                                                // intuitive interface.
    bool   fanOut     = atts.GetFanOutPoints();
    bool   doTris     = atts.GetDoTriangles();
    std::string var   = atts.GetVariable();

    if (maxSubdivs <= 1)
    {
        if (!haveIssuedWarning)
        {
            avtCallback::IssueWarning(
                "The maximum number of subdivisions must be positive");
            haveIssuedWarning = true;
        }
        return in_dr;
    }

    //
    // If we have 3D data, then return -- we only know how to execute on
    // polygonal data.
    //
    if (in_ds->GetDataObjectType() != VTK_POLY_DATA)
    {
        if (!haveIssuedWarning)
        {
            avtCallback::IssueWarning(
                 "The subdivide quadrilaterals operator only operates on "
                 "polygonal data.  You must convert this data before "
                 "using this operator.  This can be done using the "
                 "ExternalSurface operator.");
            haveIssuedWarning = true;
        }
        return in_dr;
    }
    vtkPolyData *input = (vtkPolyData *) in_ds;

    //
    // Determine what variable we are working with.
    //
    vtkDataArray *arr = NULL;
    if (var == "default")
        arr = input->GetPointData()->GetScalars();
    else
        arr = input->GetPointData()->GetArray(var.c_str());
    bool isSuitable = true;
    if (arr == NULL || arr->GetNumberOfComponents() != 1 ||
        (arr->GetDataType() != VTK_FLOAT && arr->GetDataType() != VTK_DOUBLE))
        isSuitable = false;
    if (!isSuitable)
    {
        if (!haveIssuedWarning)
        {
            avtCallback::IssueWarning(
                 "The variable to be used either was zonal or was not a "
                 " scalar or was not floating point data.  The SubdivideQuads "
                 "operator is only for nodal, scalar floating point "
                 "quantities.");
            haveIssuedWarning = true;
        }
        return in_dr;
    }

    //
    // Calculate how many new pts/cells will be needed so we can allocate
    // the correct amount of memory.
    //
    int *numQuads = new int[maxSubdivs+1];
    int *numNewPts = new int[maxSubdivs+1];
    int base = 1;
    for (int i = 0 ; i <= maxSubdivs ; i++)
    {
        numQuads[i] = base*base;
        numNewPts[i] = (base+1)*(base+1) - 4;
        base *= 2;
    }

    //
    // Make a first pass through the data to determine how much memory will
    // be needed.
    //
    vtkIdType ncells = in_ds->GetNumberOfCells();
    vtkIdType num_new_pts = 0;
    int mem_size = 0;
    vtkIdType cell_cnt = 0;
    for (vtkIdType i = 0 ; i < ncells ; ++i)
    {
        vtkIdType ct = input->GetCellType(i);
        vtkCell *cell = input->GetCell(i);
        vtkIdList *ids = cell->GetPointIds();

        bool cameFromTri = false;
        bool shouldExamine = false;
        vtkIdType id[4];

        if (ct == VTK_TRIANGLE && doTris)
        {
            //
            // We are opting to subdivide triangles.  Treat them like
            // degenerate quads because it makes everything so easy...
            //
            shouldExamine = true;
            cameFromTri = true;
            id[0] = ids->GetId(0);
            id[1] = ids->GetId(1);
            id[2] = ids->GetId(2);
            id[3] = ids->GetId(0);
        }
        else if (ct == VTK_QUAD || ct == VTK_PIXEL)
        {
            shouldExamine = true;
            id[0] = ids->GetId(0);
            id[1] = ids->GetId(1);
            id[2] = ids->GetId(2);
            id[3] = ids->GetId(3);
        }
        else
        {
            mem_size += cell->GetNumberOfPoints() + 1;
            cell_cnt++;
        }

        if (shouldExamine)
        {
            int numSubdivs = DetermineSubdivisions(arr, id, maxSubdivs, 
                                                   threshold);
            if (numSubdivs == 0 && cameFromTri)
                mem_size += 1+3;
            else
                mem_size += 5*numQuads[numSubdivs];
            num_new_pts += numNewPts[numSubdivs];
            cell_cnt += numQuads[numSubdivs];
            if (numSubdivs > 0 && cameFromTri)
                num_new_pts++; // extra pt for our degenerate quad trick.
        }
    }

    //
    // Now that we know how much memory to use, allocate all of our output
    // constructs (connectivity, points, point data, cell data).
    //
    vtkPolyData *output = vtkPolyData::New();
    output->GetFieldData()->ShallowCopy(input->GetFieldData());

    vtkPointData *oldPD = input->GetPointData();
    vtkPointData *newPD = output->GetPointData();
    vtkCellData *oldCD = input->GetCellData();
    vtkCellData *newCD = output->GetCellData();

    vtkIdType noldpts = input->GetNumberOfPoints();
    vtkIdType nnewpts = noldpts + num_new_pts;

    int ptsDataType = input->GetPoints()->GetDataType();
    vtkPoints *new_pts = vtkPoints::New(ptsDataType);
    new_pts->SetNumberOfPoints(nnewpts);
    if (ptsDataType == VTK_FLOAT)
    {
        float *old_pts_ptr = (float *) input->GetPoints()->GetVoidPointer(0);
        float *new_pts_ptr = (float *) new_pts->GetVoidPointer(0);
        memcpy(new_pts_ptr, old_pts_ptr, 3*noldpts*sizeof(float));
    }
    else if (ptsDataType == VTK_DOUBLE)
    {
        double *old_pts_ptr = (double *) input->GetPoints()->GetVoidPointer(0);
        double *new_pts_ptr = (double *) new_pts->GetVoidPointer(0);
        memcpy(new_pts_ptr, old_pts_ptr, 3*noldpts*sizeof(double));
    }

    newPD->CopyAllocate(oldPD, nnewpts);
    for (vtkIdType i = 0 ; i < noldpts ; i++)
    {
        newPD->CopyData(oldPD, i, i);
    }

    output->Allocate(mem_size);
    newCD->CopyAllocate(oldCD, cell_cnt);

    //
    // We are finally ready to do the subdivision.
    //
    vtkIdType *sub_ids = new vtkIdType[numNewPts[maxSubdivs]+4];
    vtkIdType curCell = 0;
    vtkIdType curPt = noldpts;
    for (vtkIdType i = 0 ; i < ncells ; i++)
    {
        //
        // Get the cell and see if we can subdivide it.
        //
        vtkCell *cell = input->GetCell(i);
        int ct = cell->GetCellType();
        vtkIdList *ids = cell->GetPointIds();
        bool shouldDoCell = false;
        if (ct == VTK_QUAD || ct == VTK_PIXEL)
            shouldDoCell = true;
        if (ct == VTK_TRIANGLE && doTris)
            shouldDoCell = true;
        if (!shouldDoCell)
        {
            output->InsertNextCell(ct, ids);
            newCD->CopyData(oldCD, i, curCell);
            curCell++;
            continue;
        }

        //
        // If we are going to fan out the points along the outer edge
        // of the cell, we should determine the middle, so we can use
        // it repeatedly later.
        //
        double middle[3];
        if (fanOut)
            vtkVisItUtility::GetCellCenter(cell, middle);

        //
        // We want to treat every cell like a quad for simplicity.  If it
        // is a pixel or triangle, manipulate it to look like a quad.
        //
        vtkIdType id[4];
        id[0] = ids->GetId(0);
        id[1] = ids->GetId(1);
        id[2] = ids->GetId(2);
        if (ct == VTK_TRIANGLE)
        {
            // Just double up pt 0 for now.
            id[3] = id[0];
        }
        else
            id[3] = ids->GetId(3);
        if (ct == VTK_PIXEL)
        {
            vtkIdType tmp = id[2];
            id[2] = id[3];
            id[3] = tmp;
        }

        //
        // Determine how many subdivisions to do.
        //
        int numSubdivs = DetermineSubdivisions(arr, id, maxSubdivs, threshold);

        //
        // Do some special handling for the triangle case.  Either decide
        // that we aren't doing subdivision, so create the original triangle,
        // or decide that we are and we need to create a new point for the
        // degenerate quad the triangle is going to become.
        //
        if (ct == VTK_TRIANGLE)
        {
            if (numSubdivs == 0)
            {
                //
                // Just go ahead and insert the original cell.  Since we are
                // not subdividing, there is no point in treating this triangle
                // like a degenerate quad any longer.
                //
                output->InsertNextCell(ct, ids);
                newCD->CopyData(oldCD, i, curCell);
                curCell++;
                continue;
            }
            else
            {
                // Make a new pt that is along edge id[0]-id[1].
                if (ptsDataType == VTK_FLOAT)
                {
                    float pt[3];
                    CreateEdgePoint(id, pt, 
                        (float*)new_pts->GetData()->GetVoidPointer(0));
                    new_pts->SetPoint(curPt, pt);
                }
                else if (ptsDataType == VTK_DOUBLE)
                {
                    double pt[3];
                    CreateEdgePoint(id, pt, 
                        (double*)new_pts->GetData()->GetVoidPointer(0));
                    new_pts->SetPoint(curPt, pt);
                }
                newPD->InterpolateEdge(oldPD, curPt, id[0], id[1], 0.5);
                curPt++;
            }
        }

        // 
        // Make our indexing scheme work.
        //
        int base = 1;
        int step = 1;
        for (int j = 0 ; j < maxSubdivs ; j++)
            if (j < numSubdivs)
                base *= 2;
            else
                step *= 2;

        double ratio = (double) base*step;
        int numRows = base*step+1;
        int numColumns = base*step+1;
        sub_ids[0] = id[0];
        sub_ids[numRows-1] = id[1];
        sub_ids[(numRows-1)*numColumns] = id[3];
        sub_ids[numRows*numColumns-1] = id[2];

        //
        // Set up pts inside the quad.
        //
        int I, J;
        vtkIdList *new_ids = vtkIdList::New();
        new_ids->SetNumberOfIds(4);
        new_ids->SetId(0, id[0]);
        new_ids->SetId(1, id[1]);
        new_ids->SetId(2, id[2]);
        new_ids->SetId(3, id[3]);
        for (I = 0 ; I < numRows ; I += step)
            for (J = 0 ; J < numColumns ; J += step)
            {
                bool iSpecial = (I == 0 || I == numRows-1);
                bool jSpecial = (J == 0 || J == numColumns-1);
                if (!iSpecial || !jSpecial)
                {
                    double w[4];
                    w[0] = ((ratio-I) / ratio) * ((ratio-J) / ratio);
                    w[1] = ((I) / ratio) * ((ratio-J) / ratio);
                    w[3] = ((ratio-I) / ratio) * ((J) / ratio);
                    w[2] = ((I) / ratio) * ((J) / ratio);

                    if (ptsDataType == VTK_FLOAT)
                    {
                        float pt[3];
                        CreateMidPoint(fanOut, iSpecial, jSpecial, id, w,
                            middle, pt, (float*)new_pts->GetVoidPointer(0));
                        new_pts->SetPoint(curPt, pt);
                    }
                    else if (ptsDataType == VTK_DOUBLE)
                    {
                        double pt[3];
                        CreateMidPoint(fanOut, iSpecial, jSpecial, id, w,
                            middle, pt, (double*)new_pts->GetVoidPointer(0));
                        new_pts->SetPoint(curPt, pt);
                    }
                    newPD->InterpolatePoint(oldPD, curPt, new_ids, w);
                    sub_ids[J*numRows+I] = curPt;
                    curPt++;
                }

                if (I != 0 && J != 0)
                {
                    newCD->CopyData(oldCD, i, curCell);
                    curCell++;
                    vtkIdType quad_ids[4];
                    quad_ids[0] = sub_ids[(J-step)*numRows+(I-step)];
                    quad_ids[1] = sub_ids[(J-step)*numRows+I];
                    quad_ids[2] = sub_ids[J*numRows+I];
                    quad_ids[3] = sub_ids[J*numRows+(I-step)];
                    output->InsertNextCell(VTK_QUAD, 4, quad_ids);
                }
            }

        new_ids->Delete();
    }

    output->SetPoints(new_pts);
    output->Squeeze();
    new_pts->Delete();

    avtDataRepresentation *out_dr = new avtDataRepresentation(output,
        in_dr->GetDomain(), in_dr->GetLabel());

    output->Delete();

    delete [] numQuads;
    delete [] numNewPts;
    delete [] sub_ids;

    return out_dr;
}

